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Abstract 

.^^ The aim of this paper is to study differential and spectral properties of the infinitesimal 

r \ operator of two dimensional Markov processes with diffusion and discrete components. The 

^^ infinitesimal operator is now a second-order differential operator with matrix-valued coefficients, 

T^H from which we can derive backward and forward equations, a spectral representation of the 

d probability density, study recurrence of the process and the corresponding invariant distribution. 

C All these results are applied to an example coming from group representation theory which can 

1—^ be viewed as a variant of the Wright-Fisher model involving only mutation effects. 

^ 1 Introduction 

en 

t — The connection between stochastic processes and orthogonal polynomials is very well known. For 

^ a very good account see [I^. Among them we can find diffusion processes where the state space 

C^^ is an interval of the real line, i.e. 5 C M, and the parameter set is T = [0, oo), i.e. continuous 

^Zi time. In this case, the corresponding infinitesimal operator A comes in terms of a second-order 

T-H differential operator. It is possible to calculate an expression of the transition probability density 

^ by spectral methods, i.e. in terms of the eigenfunctions of the infinitesimal operator A (see Section 

8, Chapter V of |3] or Section 13, Chapter 15 of [l5]), among other properties. In some cases, the 
invariant distribution can also be given in terms of the associated orthogonality measure of the 
C^ eigenfunctions. Prominent examples connected with spectral methods are the Orstein-Uhlenbeck 

process, population growth models or Wright-Fisher models. 

Continuous-time bivariate Markov processes with diffusion and discrete components are a nat- 
ural extension, where now the state space is two dimensional, i.e. S = S x {1,2, . . . ,N}, where 
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5 C M and A^ is a positive integer. The first component is continuous, which we will call the diffu- 
sion component, while the second one is discrete and which we will call the phase. We will denote 
these bivariate Markov processes by {Xt, Yt),t > 0. At first it would seem that this hybrid process 
is just blending a regular scalar diffusion process and a continuous time Markov process, but there 
are more subtleties specially in the case when both components are dependent each other. 

The ideas behind these coupled stochastic processes may be traced back to what is known in 
the literature as random evolutions. These were introduced by R. Griego and R. Hersh in 1969 
[7], although the seminal work was given by M. Kac several years before (see [H] for an historical 
account). The main idea of a random evolution is a model for a dynamical system whose equation 
of state is subject to random variation. Other authors like G. Papanicolaou, M. Pinsky, T. Kurtz 
or R. Varadhan, to mention a few, were involved in this field during the 1960s and the 1970s. For 
a very good account see the Hersh survey [13] and the books [l7l[l8]. In the 1980s, the Kiev school 
led by Korolyuk and Swishchuk put these processes in the framework of semi-Markov processes 
(see [5nj for a recent book on this subject). 

When the equation of state is governed by diffusion processes, we are in the case of bivariate 
Markov processes, the objects we will be interested in this paper. Now the infinitesimal operator 
A comes in terms of a second-order differential operator with matrix- valued coefficients. The basic 
definitions and properties of these bivariate Markov processes were introduced in [2j (see also Section 
p] below) . Here the author focused on the calculation of the first-passage time distribution out of 
the phases using the resolvent of the corresponding infinitesimal operator. This model was shown 
to be applicable to describing the progression of HIV in an infected individual. These processes 
have received other names in the literature, like diffusions in random environments or switching 
diffusion models (see for instance [6]). 

It is not difficult to foresee the applications that these bivariate Markov processes may have, as 
in random evolutions or switching diffusion models, in fields like population dynamics, insurance 
or financial mathematics. 

In this paper we will focus on differential and spectral properties of the infinitesimal operator 
generated by a bivariate Markov process with diffusion and discrete components. Surprisingly 
enough, this has not been explored in depth so far. The aim of this paper is to fill this gap using 
spectral analysis of differential operators with matrix-valued coefficients. 

In the first part, after some preliminary considerations in Section^ (already introduced in [2]), 
we will derive, in Section [3| the corresponding backward and forward equations, a representation 
of the (matrix- valued) probability density in terms of the eigenfunctions of the infinitesimal oper- 
ator A and some differential equations associated with some functionals, including the invariant 
distribution of these processes. 

In the second part, in Section El we will study an example (coming from group representation 
theory [16]) of these bivariate Markov processes which can be viewed as a variant of the Wright- 
Fisher model involving only mutation effects. We will use the results in Sections [2] and [3] to study 
the probabilistic aspects of this example, specially the spectral representation of the probability 
density and an explicit expression of the invariant distribution. 



2 Preliminaries 

Let {Xt,Yt),t > be a time-homogeneous Markov process assuming values in S* x {1, 2, . . . ,N}, 
where S" C M and iV is a positive integer. We use preferably boldface capital letters to denote 
matrices, boldface lowercase letters to denote vectors and standard font for scalars. We define the 
transition probability density as a matrix- valued function P{t;x,A) = (Pij {t; x , A)) , defined for 
every t > 0,x G S", and real Borel set A, as a matrix whose entry (i,j) is given by 



(2.1) 



P,,(t; X, A) = Pr{Xt e A,Yt = j\Xo = x, Yo = i}. 



For simplicity, we will assume that every entry Pij{t; x, dy) has a density, also called Pij{t; x, y), in 
such a way that we can write 



Pij{t;x,A)= / Pij{t;x,y)dy. 
J A 

Observe that P{t] x, y) is not a probability density in the common sense, but we have that all 
entries are positive and for every real Borel set A C S 



/ P{t;x,y)dyj e^ < sn, 



where ejv is iV-dimensional column vector with all entries equal to 1. The equality holds when 
A = S. 

It is assumed that P{t; x, y) satisfies the Chapman-Kolmogorov equation 



f2.2) 



P{s + t;x,y) = / P{s;x,z)P{t;z,y)dz. 
Js 



In a similar way as the classical theory of diffusion processes we have to make the following 
assumptions: 



1. For every e > 
(2.3) 



lim — P(/i; X, {y : \y — x\ > e}) = 0. 
h-s>0 h 



2. There exist N x N matrix-valued functions A{x) (with positive entries) and B{x), x G S, 
such that for any e > 



(2.4) 
and 

(2.5) 



lim — 

h^O h 



lim — 

/i-5>o h 



{y-\y-x\<e} 



{y-x)P{h;x,y)dy 



B{x) 



{y-\y-x\<e} 



{y - xfP{h;x,y)dy 



A{x). 



3. There exists an N x N matrix-valued function Q{x) = {Qij{x)) with Qij real-valued such 
that 



(2.6) 



\iin\[P{h-x,S)-I\ =Q{x). 
h^O h 



Observe that (2.3) gives continuous sample paths for the process Xt while (2.4) and (2.5) are 



matrix-valued versions of the corresponding hypothesis employed in diffusion processes (see, for 



example, [3l[T5]). Note that (2.6) is a more general form of a corresponding hypothesis employed 
in continuous time Markov chains. The formulation here is different because Q depends on x. In 
particular we have that, for fixed x, 

Qij{x) < for i = j, 
Qij{x) > for i^j; 

and for each i 

i.e., every off diagonal entry of Q{x) is a positive function and the sum of every row of Q{x) is 0, 
for every x, i.e., Q{x)e]\f = 0. 

We will denote by ^{S^) the set of all column vector- valued functions where all components 
are real-valued, bounded and Borel measurable functions on S. 
For any f{x) G 53(5"^) we define the transition operator as 



(2.7) 



Ttf{x) = E[f{Xt)\Xo 



Pit;x,y)f{y)dy. 



The transition operator (2.7) can be applied to matrix-valued functions F{x) G *B(S' ), 
with similar properties as *B(S^). We just have to perform the operator column by column. For 
simplicity we will use vector- valued analysis, but sometimes it will be more convenient to use a 
matrix-valued approach. Both are vector spaces with a norm over M. 

The family of transition operators {Tt : t > 0} has the semigroup property 



f-s+t 



T^Tf 



sJ-ti 



as a consequence of the Chapman-Kolmogorov equation (2.2). In particular, this implies that the 



transition operators commute. Therefore, the behavior of Tt near i = completely determines the 
semigroup, as in the scalar situation. 

The infinitesimal generator A of the process is defined for each smooth vector- valued function 
f{x) as 



(2.8) {Af){x) = lim I [Ttfix) - fix)] = lim ^ 



P{h;x,y)f{y)dy - fix) 



The following two important results were proved in [2]: 



Lemma 2.1 (Lemma 2.1 of [2j). A(x) and B{x), x €z S, are diagonal matrices. 

Theorem 2.1 (Theorem 2.1 of [2j). For any vector-valued function f{x),x G S, whose entries are 
real-valued, bounded and have continuous second derivatives, the infinitesimal operator A defined 



in (2.8) is of the form 



(2.9) {Af)ix) = lA{x)r{x) + B{x)f{x) + Q{x)f{x). 

For simplicity, since the coefficients A{x) and B{x) are diagonal matrices, we will denote the 
corresponding diagonal entries by crf{x) and tj(x), i = 1, . . . , N, respectively. 

Now we will give a stochastic representation of these processes, already given in Section 3 of 
[2]. We refer to Figure [T] for a better comprehension of how these processes work. 

Suppose the process {Xt,Yt) starts at state Xq = x and phase Iq = i- Let Xj be the diffusion 
process starting at the point x and evolving according to the infinitesimal generator 

(2.10) r^^^^&^+^^^^^Tx+^-^-^&o- 

XI is a diffusion process with diffusion and drift coefficients af{x) and Ti{x), respectively, and 
a killing coefficient —Qii{x) > 0. Since the diffusion XI has a lifetime tj (chosen according to 
exponential holding time of —Qii{Xt), see pp. 314 of [E]), the process is defined on [0,tj] as 
Xt = XI, Yf = i (if ti = c«, then Xt = XI and Yt = i for all i > 0). At time t = ti,Yt changes from 
phase i to another phase j ^ i. If —Q^^{Xt^) > 0, there is a probability distribution over the phases 
j 7^ i, defined as —Qij{Xt^)/Qii{Xt^),j ^ i, and Yt is selected in accordance with this distribution. 
At the moment ti, the diffusion component Xt is changed from XI to Xl, starting at XI, and 



evolving according with the generator (2.10) with index j in the place of i. The process is killed 
at the time ti + tj, where tj is the lifetime of Xl and the process is defined on [ti,ti + tj\ as 
Xt = Xl,Yt = j. If —Qjj{Xt^+tj) > 0, there is a probability distribution over the phases k ^ j, 
defined as —Qjk{Xt^+tj)/Qjj{^ti+tj)-, k ^ j, and Yt is selected in accordance with this distribution. 
Then the diffusion changes from Xf. to Xj^ at time ti + tj, where X^ starts at Xt-+t- and has 



generator (2.10) with k in the place of i, and so on. 

Observe that the sample path of the process is defined as a piecewise function and is always 
continuous for i > 0. The only difference with respect to a regular diffusion process is that this 
process can move through N different phases with different diffusion and drift coefficients in general. 
The transitions through all the phases and the lifetime spent at each phase, which depend on the 
position Xt, are derived in terms of the coefficient Q{x). 

As an illustration, in Figure [l] we display a sample path of one example of this kind of processes. 
The example is conveniently chosen so that it points out the effect of the diffusion coefficients. 

We take A^ = 3 phases and S = M. The simulation of these processes are made according to 
the stochastic differential equation 

dXt = TY, {Xt)dt + ay, {Xt)dBt, Yt = 1, 2, 3, 



where 



<rhx) 



Tiix) 



-IX, 



1,2,3. 



Here Bt denotes the standard Brownian motion. 

The three phases evolve as Orstein-Ulenbeck processes with different parameters. The diffusion 
coefficients grow with the number of phase and the drift coefficients are more intense with the 
number of phase. The transition of phases goes according to the matrix 



Q{x) 






^ ^ 400 ^ 
^2 ^ 2^ 



6(l+x2) 

_9 _ ^^ 
100, 



6(l+x2) 



1 + 



3a; 
400 



In Figure [T] we clearly see the evolution of the particle at each phase, starting at x 
2. We observe that the intensity of the process is higher with the number of phase. 



and phase 



Bivariate Orstein-Uhlenbeck process 




0.5 1 1.5 2 2.5 3 

Figure 1: Sample path of a bivariate Markov process with 3 phases. 



3 Differential and spectral properties of A 

In this section we will show the main theoretical results of this paper. We will start deriving 
the corresponding matrix-valued backward and forward equations associated with these processes. 
Then we will give an explicit formula of the matrix-valued transition density P{t; x, y) defined in 



(2.1 ) in terms of the matrix- valued eigenfunctions of the infinitesimal generator A defined in (2.9). 



Finally, we will study some differential equations associated with some functionals related with the 
recurrence and the invariant distribution of these processes. 



3.1 The matrix- valued backward and forward equations 



For all f{x) G 55(5 ) that are twice continuous differentiable we have that the matrix-valued 



expectation defined in (2.7) satisfies the backward diflterential equation 

d 



(3.1) 



dt 



Ttf{x) = ATtf{x), 



with initial condition T^f^x) = f{x), where A is defined in (2.9). This is easy to see by definition 



of A in (2.8) applied to / = Ttf, using the Chapman-Kolmogorov equation (2.2). The initial 



condition is a consequence of the assumption (2.6). We also have that Tt and A commute on all 
functions in the domain of A. 

As a consequence, the transition density P{t; x, y) satisfies the Kolmogorov backward equation 



(3.2) 



dP{t-x,y) 1 ^. ,a2p(t;x,y) 



A{x) 



dt 2 

applicable for i > and x, y € S. 



da 



+ B{x) 



dP{t;x,y) 
dx 



+ Q{x)P{t;x,y), 



The matrix-valued forward, evolution or Fokker-Planck equation for P{t; x, y) can be regarded 
as dual to ( |3.2[ ) and the pertinent variables are t and y (the state variable at time t rather than 
the initial value x). 

For any integrable vector-valued function g, define the operator 



T*g{y)= [ P*it;x,y)g{x)dx. 
Js 



Observe that the integration of this operator is with respect to x, rather than y, as in (2.7). This 
operator is adjoint to Tt in the sense that 

(3.3) {Ttf,g) = [ g*{x)( [ P{t;x,y)f{y)dy)dx= [ f g*{x)P{t;x,y)f{y)dxdy 

Js \Js / JsJs 

= j^ (j^P*{t-x,y)g{x)dx\ f{y)dy = {f^T^g). 

Here {f,g) = jag*{x)f{x)dx is a real- valued inner product for any two integrable vector- valued 
functions. Observe that in case we are using matrix-valued functions (as we will do in Section 3.2) 



we have to use the same inner product {F,G) = JgG*{x)F{x)dx, but now this inner product is 
matrix-valued. 

If / is twice continuously differentiable and vanishes outside a compact subset of S, then one may 



differentiate with respect to t in (3.3) and interchange the orders of integration and differentiation 



to get, using (|3.1|), 
(3.4) 



dt' 



/> ii^a 



d_ 
dt' 



Ttf,g) = {TtAf,g) = {Af,T*g) = {f,A*T*g). 



Here the formal adjoint A* of A with respect to the inner product (•, •) is given by 
(3.5) {A*g){y) = l^{A*{y)g{y)) - ^{B*{y)g{y)) + Q*{y)g{y). 



From (3.4) we get the forward differential equation 

d 



^T;g{y)=A*T:g{y), 



or, equivalently, for the transition density P(t;x,y) 
dP*{t;x,y) 1 9\ d 



dt 2 dy"^ ' ' dy 

In other words, taking adjoint 



A*{y)P*it-x,y)) - —{B*{y)P*{t;x,y)) + Q*{y)P*it;x,y). 



dP(t:x,y) 1 (9^ ,„, ^ ^ xx d 



(3.6) 'r = 7i7r2iPit-':-y)Mv)) - i:{P{t:^,y)B{y)) + P(tii.y)Q(y)- 



Observe that the forward differential equation (3.6) has the coefficients of the infinitesimal 
operator multiplied on the right, while the backward differential equation (3.2) on the left. 

In principle we can try to solve explicitly the backward (or forward) differential equation 
with specific boundary values to obtain a general expression of the transition probability density 
P{t;x,y). Nevertheless, this task is not straightforward at all. 

There is a very simple case, when A{x) = a'^{x)I, B{x) = t{x)I and Q{x) = Q, i.e., a constant 
matrix. In this case it is easy to see that 

P{t;x,y)=p{t;x,y)e^^, 

where p{t; x,y) is the transition probability density of the one dimensional diffusion process with 
drift t{x) and diffusion coefficient a'^{x). In this case both the continuous and discrete components 
are independent and the transition of phases and the lifetime spent on them is independent of the 
position of Xj. 

We can try to generalize this separation of variables in the same case as before, but with Q{x) 
depending on x. In this case substituting P{t;x,y) = p{t]x,y)e^^^^' in (3.2) leads to a nonlinear 
second-order differential equation satisfied by Q(x), which is not straightforward to solve. Things 
get more complicated assuming A{x) and B{x) diagonal with different entries. Even when Q{x) 
is constant it is not easy to solve (3.2). 

Therefore, we need different approaches of obtaining or approximating the transition probability 
density P{t; x, y). One of them, as we will see in the next section, is via the eigenfunctions of the 
infinitesimal operator A. 



3.2 The spectral representation of the matrix-valued transition density 

As we have seen in the last section, the explicit computation of P{t; x, y) can be something difficult 
to handle. In this section we show one way of obtaining a formula for P{t] x, y) in terms of the 
matrix-valued eigenfunctions of the infinitesimal operator A. 

Consider now the space L1^{S,C ) of all matrix-valued functions F{x) such that 



{F,F)w= f F*{x)W{x)F{x)dx < oo, 
Js 



where dW is a weight matrix with a smooth density W with respect to the Lebesgue measure (see 
[5] for a more concise definition of a weight matrix). In the above definition we mean that the 
integral is finite entry by entry. This induces a matrix-valued inner product for any two matrix- 
valued functions F, G G L^(5,C^^^), denoted by 

(3.7) {F,G)w= [ G*ix)Wix)F{x)dx. 

Js 

This is not an inner product in the common sense, but it has properties similar to the usual scalar 
inner products. It is also possible to define a scalar product between two matrix- valued functions 
(see [1]), given by 

{F,G) = Tr{{F,G)w)- 

Therefore, L'^{S,C^''^) with the norm \\F\\w = Tr ((F, F)h^)^/^ is a Hilbert space and (3.7) is 
the inner product (in fact, the set of equivalence classes F ^ G ii \\F — G\\w = 0). 

The idea behind this method is to find a representation of the transition operator TtF[x) (see 



(2.7)) using the method of separation of variables as a superposition (linear combination) involving 
eigenfunctions of the infinitesimal operator A. 

Suppose that we know a set of countable normalized matrix-valued eigenfunctions ($„(x))„ of 
A with corresponding eigenvalues Vn (Hermitian), i.e. A^n{x) = ^n{x)Tn- Normalization means 
{*^n-, *m)w = ^nml for Certain weight matrix W . 

If the set of finite linear combinations of eigenfunctions is complete in ^^(5, C^^^), then each 
F G L'y^{S, C ) has a Fourier expansion of the form 

oo 
F{x) = Y,^n{x){F,^n)w■ 
Consider the superposition defined by 

oo 

UF{t,x) = ^*„(x)e*i^"(F,*„)H^. 

n=0 

Here the exponential for matrix- valued argument is defined using the usual Taylor expansion. It is 
easy to see that UF{t,x) satisfies the backward differential equation 

d 

— UF{t,x) = AUFit,x), 



with initial condition Uf{0,x) = F{x). So if there is uniqueness for a sufficiently large class of 
initial functions F then we get 

TtF{x) = UF{t,x), 

which means 

I P{t;x,y)Fiy)dy = J f f] *„(x)e*r"$;(y)l^(y) J F{y)dy. 

Therefore, P(t; x, y) is given by 

oo 

(3.8) P{t- x,y) = Y, ^n{x)e'''-K{y)W{y). 

n=0 

Remark 3.1. It is well known that ^n(x) = ^n{x)Un with U*Un = I is always another family 
of normalized eigenf unctions with eigenvalues U*TnUn- However, it is straightforward to see that 



the spectral representation (3.8) of P{t;x,y) is invariant under this change. 



Observe that the condition of having eigenfunctions with Herniitian eigenvalues forces the in- 



finitesimal operator A to be self-adjoint with respect to the inner product (3.7), i.e. {AF, G)w = 
{F,AG)w for all F,G € L|^^(S', C^^^). This means, using the formal definition of the adjoint of 



A given in (3.5), that the coefficients of the infinitesimal operator A are subject to certain symmetry 



equations with respect to the weight matrix W (with the corresponding boundary conditions) 

A*{x)W{x) = W{x)A{x), 
/g g^ B*ix)W{x) = {W{x)A{x)y - W{x)B{x), 

Q*{x)W{x) = hw{x)A{x))" - {W{x)B{x)y + W{x)Q{x). 

These symmetry equations were already derived in [5^ ilO]. The goal in those papers were to 
find matrix- valued orthogonal polynomials satisfying second-order differential equations of Sturm- 
Liouville type (i.e. the differential coefficients are matrix polynomials of degree less than or equal 
to the order of differentiation). Observe that our infinitesimal operator A does not have to be in 
principle of Sturm-Liouville type. In fact, we are more interested in examples with Q{x) depending 
on X, along with diagonal coefficients A(x) and B(x). In this line, we are closer to differential 
operators as the ones introduced in [9] . Later it was shown (see [11] ) that the operators introduced 
in [9] are intimately related to second-order differential equations of Sturm-Liouville type. 

3.3 Differential equations associated with some functionals 

In this section we will study some differential equations associated with some functionals related 
with the recurrence and the invariant distribution of these processes. 

The hitting time of a one dimensional diffusion process {Xt, < t < C} to the state z is defined 
by Tz = inf{t > : Xt = z} or T^ = oo ii Xt ^ z for < t < (. We use the notation 

T* = min{ TcTrf}, 



10 



for the hitting time to c or d, the first time Xt = c or Xt = d, where c,d£ S = {a,b). For processes 
starting at Xq = x G (c, d), this is the same as the exit time of the interval (c, d). We have to 
accommodate this definition with our bivariate Markov processes. The solution will be introducing 
matrix-valued analysis. 

We will assume that our process is regular in the interior of S, i.e. 

PrjTj, < oo, Yt^ = j\Xo = x,Yo = i}> 0, 

for any x,y in the interior of S and any phases i,j. 

Consider the matrix- valued function U{x) = {Uij{x)), where every entry (i,j) is given by 

(3.10) Uij{x) = PrjTrf < r„ Yt, = j\Xo = x,Yo = i}, c < x < d, 

i.e., starting at a; G (c, d) and phase i, Uij{x) is the probability that the process reaches d before c 
and at that time the process is in phase j. 



Theorem 3.1. The matrix-valued function U{x) defined in (3.10) possesses two bounded deriva- 



tives for c < X < d and satisfies the following differential equation 

(3.11) ^A{x)U"{x) + B{x)U'{x) + Q{x)U{x) = 0, c < x < d, U{c) = 0, U{d) = I. 

Proof. The proof is essentially the same as in the scalar case. For a heuristic justification see pp. 
193 of |15j . For more technical details the reader should consult Proposition 9.1 of Chapter V in 

m- □ 

Remark 3.2. In |3 \15\j we can find an explicit solution of this problem in the scalar situation 
(where Q{x) = 0) in terms of what they call scale function and speed density of the process. The 
presence of Q[x), not zero in general, leads to a more complicate system of second- order differential 
equations with initial conditions 

1 ^ 

(3.12) -ct2(x)C/^^(x) + Ti{x)U[^{x) + Q^^{x)Uij{x) + J^ Q^k{x)Ukj{x) = 0, i, j = 1, . . . , TV. 

With the above considerations we can define recurrence and transience of bivariate Markov 
processes (in the same line as Definition 9.1 of Chapter V in [3j). Write the matrix- valued function 
Rx,y = {Rx,y)i,j where every entry (i, j) is given by 

{Rx,y)i,j = Pr{Xt = y,Yt = j\Xo = x,Yo = i}, x,y £ S, 

i.e., the probability that starting at x and phase i, the process ever reaches y and phase j. 



3.11 ) letting c ^>- a and d ^>- y 
letting c ^ b and d ^ y. 



For y > X, Rx,y will be the solution of the differential equation (t 
and, for y < x, Rx,y will be the solution of the differential equation (3.11 ) 

Definition 3.1. A state y is called recurrent if Rx^ye]\f = e^r for all x G S, where e^ is N- 
dimensional column vector with all entries equal to 1, and transient otherwise. If all states in S 
are recurrent, then the bivariate Markov process is said to be recurrent. 



11 



Consider now, for a continuous matrix- valued function G{x), finding V{x) given by 

rT* 



(3.13) 



V{x) = E 



I G{Xs)ds\XQ = X 
Jo 



c < X < d, 



where E is defined in (2.7). 



Theorem 3.2. The matrix-valued function V{x) defined in (3.13) possesses two bounded deriva- 
tives for c < X < d and satisfies the following differential equation 

(3.14) ^A{x)V"{x) + B{x)V'{x) + Q{x)V{x) + G{x) = 0, c < x < d, V{c) = V{d) = 0. 

Proof. As before, the proof is essentially the same as in the scalar case. For a heuristic justification 
see pp. 194 of [15j. For more technical details the reader should consult Proposition 10.1 of Chapter 

V in 131. D 



In particular, in the case where 



/I 1 

1 1 



1\ 
1 



(3.15) G{x) = ENel, = 

V 1 ••• 1/ 

we have that the matrix- valued function V{x) = {Vij{x)) gives 

(3.16) Vij{x) = E [T\Yt* = j|Xo = x,YQ = i], 

i.e., starting at Xq = x and phase i, Vij{x) is the mean time to reach either c or d and at that time 
the process is at phase j. 

Remark 3.3. Again, in PJ.IJ5]/ we can find an explicit solution of this problem in the scalar situation 
in terms of the Green function of the process, but the presence of Q{x) will lead to a system of 



differential equations similar to (3.12), which is more difficult to solve 



Definition 3.2. Assume our process is recurrent. By letting c ^ a and d ^>- b in (3.14), we say 



that the bivariate Markov process is positive recurrent if all entries ofV{x) in (3.16), with G{x 



given by (3.15), are finite. Otherwise is null recurrent. 



Finally, let us discuss the invariant distribution of a bivariate Markov process, which is the 
behavior of P{t;x,y) as t — )■ oo. This distribution should be independent of the initial state 
Xq = X and the initial phase Iq = i- Therefore we should expect a row vector-valued invariant 
distribution 

tp{y) = {i^i{y),i^2{y), ■ ■ ■ ,ipN{y)), 



satisfying 



^(y) = / ^{x)P{t;x,y)dx, t > 0. 
Js 
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It is also expected, as in the scalar situation (see Theorem 12.2 of Chapter V in [3j or pp. 241 of 
|15j). that the existence and unicity of such invariant distribution is restricted to the case when the 
process is positive recurrent. 



It is easy to see, mimicking the derivation of the forward differential equation (3.6), that il^{y) 
satisfies 

(3.17) ^itP{y)A{y))" - my)B{y)y + ^(y)g(y) = 0. 

In order to have a vector- valued probability distribution we must impose that < ipj{y) < 1 and 

(3.18) (j ^l^{y)d^ CN = 1. 

If the infinitesimal operator A is self-adjoint with respect to certain weight matrix W then we 



can solve directly (3.17). This equation resembles the third symmetry equation in (3.9). Substi- 
tuting 

i,{y) = ce*^W{y), c> 0, 



in (3.17) and using the third equation in (3.9) give 
^1 



ce 



N 



-(W{y)A{y))" - {W{y)B{y)y + W{y)Q{y) 



cel,Q*iy)Wiy) = 0, 



as a consequence of Q{y)eN = 0, for all y. The normalization constant c is given by (3.18): 

1 

e*NW{y)eNdy 

Therefore an explicit formula for the invariant distribution in this case is given by 

-1 



(3.19) 



ip{y) 



I e*^W{y)eNdy] e^T^(y). 



4 A variant of the Wright-Fisher diffusion model involving only 
mutation effects 

The Wright-Fisher diffusion model involving only mutation effects considers a big population of 
constant size composed of two types A and B. The intensity of mutation from A to B and from 
B to A are given by two positive constants, -^ and ^^^, respectively (a,/3 > —1). As the size 
of the population tends to infinity, it is very well known (see pp. 177 of [15J) that this model can 
be described by a diffusion process whose state space is the unit interval S = [0, 1] with drift and 
diffusion coefficient given by 



(4.1) 



t{x) = a + 1- x{a + P + 2), a^{x) = 2x{l - x), a,/3>-l. 
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The variable x is the rate of the number of A-types with respect to the total population. These 
coefficients correspond with the Jacobi diffusion model, and the probability density can be described 
in terms of the Jacobi polynomials on [0, 1]. 

In this section we will study a variation of the Wright-Fisher model with only mutation effects 
given by a bivariate Markov process whose state space is [0, 1] x {1, 2, . . . , N}, so that we can apply 
all the analysis introduced in the previous sections. In our case, there will be A^ possible phases 
with different drift coefficients (the diffusion coefficient will always be the same). Also, there will 
be an extra positive parameter k, depending on /3, which measures how the process moves through 



all the phases. All these considerations will be studied deeply in Sections |4.2| and [4~3 

The example we study comes from group representation theory, and it is related with matrix- 
valued spherical functions associated with the complex projective space. It was first introduced in 
[9] for special choice of the parameters (see also \TT\ ) . Later in \T6\ it was related with matrix- valued 
orthogonal polynomials. In Section [4.1| we will give the coefficients relevant in this example to treat 
it as a bivariate Markov process. The recurrence coefficients of this example have also been related 
with certain quasi-birth-and-death processes in |8|. Recently, new applications of these processes 
have been related to urn or Young diagram models |12j . In both cases the differential operator 
played no role. In this paper we will give a meaning of the differential operator in terms of bivariate 
Markov processes. 

4.1 A family of examples arising from the complex projective space 

In what follows we will use Eij to denote the matrix with 1 at entry (i,j) and otherwise. Let 
N € {1,2,...}, a, (3 > —1 and < k < /3 + 1. M will denote the nilpotent upper triangular matrix 



M ='^E, 



N-l 



1=1 

while J and J the following diagonal matrices 

N N 



J = Y,iN - i)Eii, J = Y,ii - ^)Eii = {N- 1)1 - J. 



i=l i=l 

Finally, the diagonal matrix H will denote 

N 



H — y ^LOjEii, 

i=l 

where 

fp-k + i-l\ fN + k-i-1 

"^=V ^-l )[ N-i 

Observe that the weights oji correspond to the weights of the Hahn polynomials supported on 
{1, 2, . . . , N} with parameters /3 — k > —1 and k — 1 > —1 (see for instance pp. 345 of ^j). 
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The pair {"H^,^}, where 

(4.2) W{x) = x''{l-xfHx-^, a;e(0,l), 

and 

Id? d d° 

where 

(4.4) A{x) = 2x{l - x)I , B{x) = {a + l)I + J-x{{a + (3 + 2)I + J), 

Q{x) = — *— (j ({13 -k + 1)1 + j){M-I) + x (j{kl + J){M* - I 

is a symmetric pair, meaning that A is self-adjoint with respect to W (with respect to the inner 



product (3.7)). See [T0l[T6] for definitions. 



The tridiagonal matrix Q{x) can be written as 

N N N-l 



(4.5) Q{x) = ^ Hi{x)E^_i - ^{Xi{x) + fj,i{x))Eii + ^ Xi{x)Ei^i+i, 

where 



i=2 j=l j=l 



(4.6) 



Xi{x)=-^{N-i){i + p-k), 
1 — X 

fii{x) =— ^(i - 1){N -i + k). 
1 — X 



We remark that Xi{x) and /ij(x) are x-dependent variations of the coefficients of the second-order 
difference equation satisfied by the Hahn polynomials (see for instance pp. 346 of [1]). 

Observe that A{x) and B(x) are diagonal matrices and the diagonal entries of A(x) are positive 
for all X £ (0, 1). Also Xi{x) and fii{x) are always positive and Q{x)eN = for every x G (0, 1). 
Therefore Q{x) is the intensity matrix of a classical continuous time Markov chain with finite 
state space. All these requirements fit perfectly with the definition of bivariate Markov processes 
introduced in Section [21 

We have many formulas and properties associated with this example at our disposal in the 
literature. In particular we can find a set of normalized eigenfunctions ($n(x))„ of the differential 



operator A (4.3) with the corresponding eigenvalues (r„)n- This will allow us to get a formula for 



the matrix- valued transition density P{t; x, y) of this process using formula (3.8). 

The eigenfunctions ($,i(x))„ of ^ (also known as the matrix-valued spherical functions associ- 
ated with the complex projective space, see [9]) are very close related to matrix- valued orthogonal 
polynomials and they have been studied in several papers. We use the family of matrix-valued 
orthogonal polynomials given in [8], since it contains most of the structural formulas we need for 



this paper. First we describe the relation between the symmetric pair {Vl^,^} given in (4.2) and 
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(4.3) and the symmetric pair {PV^,^} given in ^. This transformation was introduced for the first 



time in [TT] (see also [16]). Denote ^(x) the following matrix 

where T is the constant upper triangular matrix given at the beginning of Section 7 in [S]. 
The relation between the pair {W^,^} and {V^,^} is the following 

(4.7) W{x) = **(x)Vr(x)*(x), AF{x) = *"^(x)^(*(x)F(x)) . 

Take now the family of matrix- valued orthogonal polynomials introduced in [8], which we will 
denote by {Rn{x))n (in [8j they are denoted by ((55^(x))„). Then the family {Rn{x))n is eigenfunc- 
tion of the differential operator A, i.e. ARn{x) = i?„(x)r„, where r„ is the diagonal matrix 

(4.8) r„ = -n^I - n{{a + /3 + N)I + J) - J{{a + /3 - A: + 1)1 + J). 



This family has very special properties. In particular it is very easy to normalize since the norms 
^||~, given in formula (7.3) ( 

Therefore, define the family 



|fin||~, given in formula (7.3) of [H], are diagonal matrices. 



^n{x) = ^{x)Rn{x)\\Rn\\^. 

We have that this family is eigenfunction of the infinitesimal generator A, i.e. 

J\'Syi\X) — VPji(Xjl^, 



with Vn defined in (4.8), and orthonormal with respect to W{x) (4.2). Observe that ($„(x))„ are 



also matrix- valued polynomials, but not of degree exactly n (more precisely of degree n + N — 1). 



Therefore, according to formula (3.8), we have 



(4.9) P{t; x,y) = Y, *n(^)e*^"*:(y) W(y), 

n=0 

or, in terms of the pair {W, A}, 

P{t;x,y) = *(x) [ f;i^„(x)e*^"||i^„||~2i^;(y)l^(y) j *-i(y). 

\n=0 / 

4.2 Probabilistic interpretation 

Consider the Wright-Fisher diffusion model involving only mutation effects for a big population of 
constant size composed of two types A and B. The intensity of mutation from A to B and from B to 
A are described at the beginning of Section [4] and are given in terms of the parameters a,/3 > —1 
of the example introduced in Section |4.1[ Now there will be an indicator denoted by the extra 
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parameter k, depending on the intensity /3, which helps the population of A to survive against the 



population of B. We will give more details in Section 4.3 



First we remark the importance of having a representation of the matrix-valued probability 



density P{t;x,y) given by (4.9), which can be very well approximated by the first eigenfunctions 
^n{x) and the corresponding eigenvalues r„. We can calculate very accurately what is the prob- 
ability of reaching any Borel set B C (0, 1) and certain phase Yt = j in time t given that initially 
we start from state Xq = x and phase Yq = i. As an illustration, for N = 4 phases, starting at 
X = 1/2, at time t = 1 we have 



Prixi G 



Xo 



/0.12410905 0.08138740 0.08920446 0.1633878\ 

0.11006872 0.07334748 0.08181737 0.1527569 

0.09668381 0.06555764 0.07453306 0.1420720 

VO.08394494 0.05801744 0.06735379 0.1313385/ 



where the entries (i,j) represent the initial phase Yo = i and final phase Yi = j. Here we have 
taken the values of the parameters a = /3 = and k = i^- In this case only 3 eigenfunctions are 
needed in order to get a very accurate approximation of the probabilities. 

Our family of bivariate Markov processes (Xj, 1^) with state space [0, 1] x {1, 2, . . . , A^} can be 
described in terms of the stochastic differential equation 

dXt = TY,{Xt)dt + a{Xt)dBt, 



where (see (4.4)) 

(4.10) Ti{x) = a + l + N -i-x{a + l3 + 2 + N -i), a'^ {x) = 2x{l - x) , Yt = i = I, . . . ,N. 



The description of how the process moves through the different phases is given by Q{x) ( |4.5[ ). 
Observe that Q{x) is tridiagonal, meaning that transitions of phases are only allowed between 
adjacent phases (like a birth-and-death process). 

The behavior of Xt depends on the phase Yt and the transition of phases and the waiting times 
at each phase Yj depend also on the position of Xt. About the continuous part Xt we are interested 
in studying the behavior at the boundaries values, in this case and 1. On the contrary, from the 
discrete part Yt, we will be interested in analyzing how the process moves through the different 
phases, including the waiting times at each phase. The behavior of the bivariate Markov processes 
{Xt,Yt) will be a combination of both aspects. 

The diffusion coefficient o"^(x) is always the same, but the drift coefficients Tj(x) depend on the 
phase i. For a given phase i, moving to the next phase i + 1 means decreasing one unity of the 
intensity of mutation B ^ A {m. terms of the parameter a), and moving to the previous phase i — 1 
means increasing one unity of the same intensity of mutation B ^ A. The limiting situation is at 
phase N, where we recover the regular Jacobi diffusion (4.1). The intensity of mutation A ^ B 



(in terms of the parameter /3) is not affected in these coefficients. 

In page 239 of |.15J one can find a classification of the boundary states and 1 of the Wright- 
Fisher model involving only mutation pressures for Jacobi diffusions in terms of the parameters. 
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Both states and 1 are regular or absorbing boundaries when — 1 < q,/3 < 0, respectively. While 
both states and 1 are entrance or reflecting boundaries when a, /3 > OJj 

For our bivariate Markov process observe that the boundary state 1 preserves this behavior for 
any phase. But the boundary state only preserves this behavior at phase N. For the rest of 
phases i G {1,...,A^ — 1} the boundary state is always reflecting since a + N — i > (see (4.10 )). 
Therefore 

is an absorbing boundary for — 1 < /3 < 0, 

is a reflecting boundary for /3 > 0. 



1 







is an absorbing boundary 



for 
for 



1 < a < 
a > 0. 



AND phase N, 



is a reflecting boundary 

In Figure [2] we can check this behavior for four different situations. In all of them we have taken 
A^ = 3 phases. In this figure we have not marked with numbers the transitions of phases (in vertical 
lines) since we are only interested in the behavior at the boundaries and 1. In the first plot the 
parameters are taken so that the boundaries and 1 are both reflecting, meaning that the sample 
path is never going to reach or 1. In the second plot, 1 is absorbing, while is reflecting. We 
observe that 1 can be reached in principle at any phase. In the third plot, on the contrary, is 
absorbing while 1 is reflecting. In this case the only phase where the sample path reaches is only 
on the last third phase. In the rest of phases, becomes reflecting. Finally, in the fourth plot, the 
parameters are taken so that and 1 are both absorbing. The process can reach boundary 1 at 
any phase, but, again, boundary only when the process is at the third phase. 



According to Definition 3.1 our process is always recurrent for a,/3 > 0, while transient other- 
wise. This fact has been checked by looking extensively at the numerical solutions of Rx,y via the 



second-order differential equation (3.11). In a similar way, we can state, using Definition 3. 2l that 



our process is always positive recurrent for a, /? > by looking the numerical solutions of V{x) via 



the second-order differential equation (3.14) with G{x) given by (3.15). 

Let us now study a couple of aspects from the discrete component Yj. First we will study the 
waiting times at each phase depending on the parameters and secondly we will study the tendency 
of moving forward or backward in phases. All these facts depend on the position of Xt. 



For the waiting times we have to look at the diagonal entries of Q{x) given by (see (4.5) and 



4.6)) 








'Qn{x) = 


-^{N-l)i/3-k + l) 


« 


Qii{x) = - 


-^[{N-i){i + P-k) + x{i 




^Qnn{x) = 


--j^AN-l)k 



l){N-i + k)] 



for i - 


= 1, 


for i - 


= 2,. 


for i - 


= A^. 



,A^-1, 



We have that if x is close to 1, the diagonal coefficients are very large, meaning that all phases 
are instantaneous, i.e., the waiting times at each phase are very short until x is far from 1. This 
can be checked in the graphs of Figure [2j If x is close to 0, we have that Qnn{x) is very small, 
meaning that phase A^ becomes absorbing, i.e., if the process enters this phase and the position of 



^A reflecting boundary cannot be reached from tlie interior of the state space, but it is possible to consider 
processes that begin there. An absorbing boundary can both be entered and left. 
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1 , |3=1 , k=1 .5, 20 changes 




01=1 , p=-0.4, k=0.3, 30 changes 




0.2 0.4 



0.6 0.8 1 



a=-0.4, |3=1 , k=1 .5, 20 changes 



a=-0.4, |3=-0.5, k=0.25, 30 changes 





Figure 2: Different behaviors of the boundaries and 1 for different values of a, /3 and k. In the 
first plot and 1 are reflecting. In the second plot, 1 is absorbing and reflecting. In the third 
plot, 1 is reflecting and is absorbing (only at phase 3). Finally in the fourth plot, (only at phase 
3) and 1 are absorbing. 



Xt is close to 0, then it tends to spend long periods of time there (see the last two plots of Figure 
^. Also phase N becomes absorbing (no matter of the position of Xt) if the parameter k is close 
to (see first plot in Figure [3]). The middle phases (from 2 to A^ — 1) are never absorbing, since 
Qii{x),i = 2, . . . ,N — 1, are never close to 0. Finally, if k is close to /3 + 1, we have that phase 1 
is absorbing (see second plot in Figure |3]). For the rest of cases, the value of these diagonal entries 
and the position of Xt determine the waiting time in that phase. 

For the tendency we have to study, given that the process is at phase Yt = i, the probability of 
moving to the next or to the previous phase. Obviously, when the process is at one of the boundary 
phases 1 or A^ it can only go forward or backward, respectively. This probability is given by 



Pr{yi 
Pr{Yt 



Yt 
Yt 



i + 1} 



Xi{x) 



Xi{x) + Hi{x) 

\i{x) + fMiix) 



where the coefficients Xi{x) and Hi{x) are given by (4.6). Here x is the state Xt of the process at 
the time t of transition of phases. The forward or backward tendency depends on if /ij(x) < Xi{x) 
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a=1 ,|3=2,k=0.2, 7 changes 




a=0.5,p=2,k=2.8, 10 changes 




Figure 3: Different behaviors of the phases 1 and 3 for different values of k {N = 3). In the first 
plot k is close to 0, so the process tends to spend more time in the last absorbing phase 3. In the 
second plot, k is close to /? + 1, so the phase 1 is absorbing. 



or fii{x) > Aj(x), respectively. Define the threshold values 
(4.11) 



(N-i)ii + (3-k) . 
xo[i) = — — TTTT^ : —, ^ = 2, 



,N-l. 



{i-l){N-i + k)' 

These values are the x— solutions of the algebraic equation Xi{x) = fii{x). Therefore, at the moment 
of changing from phase Yt = i to the next or previous phase, we will have a forward tendency when 
X < XQ{i), and a backward tendency when x > xo{i). In order to be relevant, this threshold value 
xo(i) must be less than 1. Otherwise, there is always a forward tendency. 

Let us study the cases for which the threshold values have a meaning. If xq^i) > 1 for all 
i = l,...,A^ — 1, then we will always have a forward tendency. This happens when 

(/3 + 1) 



k< 



N -1 



Nevertheless, if 
(4.12) 



k> 



(/3 + l)(iV-i) 



N -1 



for some i there can be forward and backward tendency at the same run, depending on the position 
of Xt at the moment of transition. We have then a maximum backward tendency for i = 2 in 
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/3+1 
N-1 



2(/3+l) 
N-1 



3+l){N~3) (/3+l)(Af-2) 



N-1 



N-1 



/3 + 1 



Figure 4: Range of A; for the study or the tendency of the process. 



(4.12). The range of the values of k depending on the phase Yt = i are described in Figure El 



Therefore we have 



Maximum forward tendency for 
Maximum backward tendency for k > 
Backward/forward tendency otherwise 



'*' ^ N-1 ' 

W+l){N-2) 
N-1 ■ 



This also explains why phases 1 or A^ are absorbing in terms of the values of k. 

Let us study closer this tendency behavior and the threshold values ( 4.11[ ) for an specific case. 
We fix A^ = 5 phases and the parameters a = |,/3 = 1. The interval (0,2) in Figure 4 is split 
in four subintervals (0, 1/2] U (1/2, 1] U (1, 3/2] U (3/2, 1). We study four values of k taken in the 
middle of each subinterval. In Figure [5] we have plotted the threshold function (4.11) (depending 
on the phase i) for these four values of k. We observe that the threshold value at the boundary 
phase 5 is 0, meaning that the next phase of the process is always 4, and the threshold value at the 
boundary phase 1 is oo, meaning that the next phase of the process is always 2. 




2.5 3 3.5 

Phases 



Figure 5: 
i/3 = l. 



Threshold values (4.11) (in big dots) for different values of k for N = 5 phases, a 
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We observe, for k G (0, 1/2) in Figure pi that there is no threshold between and 1, i.e. always 
a forward tendency. In the following table we can confirm this behavior for a sample path of the 
process. Here in the first row we have written the sequence of transitions of phases and in the 
second row the position of Xt at the end of that phase. 



Phases 


1 


2 


3 


2 


3 


4 


5 


4 


5 


4 


5 


4 


5 


4 


5 


Position 


.68 


.79 


.82 


.80 


.69 


.71 


.76 


.55 


.86 


.86 


.85 


.71 


.59 


.76 


.69 



For k G (1/2, 1) in Figure^ only phase 4 has a threshold value (around 0.8). That means that 
for phases 2 and 3 there is always a forward tendency. At the moment of transition at the end 
of phase 4, if the position of Xt is greater that the threshold value 0.8, then there is a backward 
tendency. Otherwise, there is a forward tendency. In the following table, as before, we have run a 
sample path for this situation: 



Phases 


1 


2 


3 


2 


3 


4 


3 


4 


5 


4 


5 


4 


5 


4 


5 


Position 


.49 


.53 


.56 


.51 


.46 


.87 


.52 


.45 


.12 


.13 


.20 


.04 


.13 


.12 


.37 



For k G (1,3/2) in Figure [sl phases 3 and 4 have a threshold value (around 0.85 and 0.55 
respectively). Again, phase 2 has a forward tendency, while phases 3 and 4 depend on the position 
of Xt as before. In the following table we have run a sample path for this situation: 



Phases 


1 


2 


3 


4 


5 


4 


3 


2 


3 


2 


3 


4 


5 


4 


3 


Position 


.71 


.83 


.78 


.43 


.83 


.81 


.88 


.78 


.64 


.86 


.69 


.72 


.42 


.58 


.31 



Finally, for k G (3/2, 1) in Figure pi all middle phases have a threshold value (around 0.8, 0.6 
and 0.4 respectively). This is the situation with the maximum backward tendency, but it depends 
strongly on the position of Xt- Again, in the following table we have run a sample path for this 
situation: 



Phases 


1 


2 


1 


2 


3 


2 


3 


4 


3 


2 


1 


2 


3 


4 


3 


Position 


.56 


.72 


.43 


.47 


.75 


.66 


.64 


.87 


.55 


.62 


.65 


.17 


.57 


.62 


.43 



Finally, we will give the explicit expression of the vector- valued invariant distribution i/'(y) (see 
the end of Section^ for definitions). We can get an explicit formula in this case since .4, given by 



4.3), is self-adjoint and we have an explicit expression of the orthogonality measure W . Formula 



3.19) gives 



(4.13) 



^(y) 



j e*^W{y)eMdy\ e*^W{y). 



where W is given by (4.2). A closer look to the 0-norm and H shows, for j = 1, . . . , A'^, 

'a + /? + A\ (/3 + N){k)N-j{l3 -k + l)j_i 



(4.14) 



V'j(y) = y' 



a+N-j 



(1-y) 



/3 



A-1 



,J-lyV " / {a + p-k + 2)N-i 

where {a)n denotes the Pochhammer symbol defined by {a)n = a{a + 1) • • • (a + n — 1) for n > 0, 



(a)o = 1. Another possibility of obtaining (4.13) is by looking the spectral representation of 



P{t;x,y) (4.9). From this representation we observe that the limit as t — )• oo only depends on the 
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eigenvalues Tn, given in (4.8). All diagonal entries of r„ are strictly negative for the range of all 



parameters n, a, /3, k and A^ involved, except for n = and entry (1, 1) of Fq, which is 0. Therefore 



taking limits in (4.9) when t — ?■ oo, we have that the only relevant element is the first term of the 
sum, given by 

lim Pit;x,y) = ^o{x)En^Uy)Wiy). 



t— >oo 



An straightforward computation of the expression above, using the notation given in Section 4.1 



leads to a matrix where all rows are equal to (4.13). 



This vector- valued invariant distribution (^4.13) is valid only when the process is positive recur- 



rent, i.e. a,/3 > 0. For — 1 < a,/3 < 0, (4.13) is also meaningful, but the boundary points of the 



process are now absorbing, meaning that the correct vector-valued invariant distribution of such 



cases involves mass jumps at the boundaries and 1 plus a density portion of the form (4.13). 



As an illustration, in Figure [gI we have plotted ipj{y), j = 1, ... A, in ( |4.14 ) for A = 2, 3, 4, 5. 



We have fixed the parameters to a,/3 = 1,A; = 5/4. We remark that Jq ipj{y)dy < 1 for each 

component, but Ylj=i /o '^j{y)dy = 1- 

From these plots we immediately see that, for a large time, if the process is at one state in 
(0, 1/2), then there is a very high probability that the process is located at one of the last phases. 
Nevertheless if the process is at one state in (1/2, 1), there is more or less the same probability that 
the process is located at any phase (higher probability for the first phases). 

This interpretation may change depending on the value of k. If k is close to /3 + 1 then we have 
the maximum backward tendency and the process is more likely to be at one of the initial phases. 
Nevertheless, if k is close to 0, then we have maximum forward tendency, and the process is more 
likely to be located at one of the final phases. 

4.3 Final remarks 



From the considerations in Section 4.2 we are ready to give an interpretation of this bivariate 
Markov process in relation with the Wright-Fisher diffusion model involving only mutation effects. 
Consider a big population of constant size composed of two types A and B, where the intensity of 
mutation from A to B is given in terms of the parameter /3 and the intensity of mutation from B to 
A is given in terms of the parameter a. The value that the process Xt takes in the interval (0,1) 
is the rate # A/(# A+# B), i.e. if Xt < 1/2 then # B> # A and if Xt > 1/2 then # A> # B. 
We assume that the intensity of mutations a and /3, are more or less the same. If a > /3 then the 
population of A's tends to survive against the population of B's (the opposite for a < /?). 

Now the process can move through N different phases. At phase A^ the process behaves as 
a regular Wright-Fisher model (in terms of the Jacobi diffusion). At phase N — 1 the intensity 
B ^ A increases in one unity, at phase N — 2 the intensity B ^ A increases in two unities, and so 
on, until the first phase, where the intensity B ^ A increases in A^ — 1 unities. Therefore moving 
through the first phases means that the intensity of mutation B ^ A is much higher if we compare 
with the intensity of mutation A ^ B, i.e. it helps to increase the population of A's against B's. 

The meaning of the indicator parameter k (depending on /3) is the following. If k is close to 0, 
then we have a forward tendency, i.e. the process tends to spend more time at the last phases, i.e. 
the populations of A's and B's "fight" more or less in the same conditions (remember that we are 
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Figure 6: The components of the vector- valued invariant distribution ipj{y),j = 1,...,N (see 
( |4l4l )), for iV = 2,3,4,5, anda,/3 = l,/c = 5/4. 



assuming that a is more or less equal to /?). On the contrary, if this indicator k is close to /3 + 1, 
then we have a backward tendency, i.e. the process tends to spend more time at the initial phases, 
i.e. it helps the population of A's to survive against the population of B's. For the middle values 
of k there can be forward/backward tendency, and it depends on the number of A's and B's, i.e. 
the position of Xt- 

5 Conclusions 

In this paper we have studied bivariate Markov processes with diffusion and discrete components 
from its matrix- valued infinitesimal operator properties. This infinitesimal operator is a second- 
order differential operator with matrix-valued coefficients with certain properties. We have given 
some theoretical results related with these processes, like the backward/forward equations, the 
spectral representation in terms of the eigenfunctions of the infinitesimal operator, recurrence 
considerations and the invariant distribution. All these results have been applied to an example 
coming from group representation theory and related with matrix- valued spherical functions, which 
can be viewed as a variant of the Wright-Fisher diffusion model involving only mutation effects. 
Since we have an explicit set of eigenfunctions with the corresponding orthogonality measure we 
can approximate very accurately the matrix-valued probability density P{t; x, y) and get an explicit 
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expression of the vector-valued invariant distribution, among other properties. 

Certainly, new examples of bivariate Markov processes may be possible, specially when the 
diffusion state space is other than a bounded interval. In [S] one can find other types of examples. 
The example studied in this paper, also known as the one-step example, is just one class of a very 
large class of matrix-valued spherical functions associated with the complex projective space. The 
second-order differential operators associated with these matrix- valued spherical functions seem to 
fit very well in the theory of bivariate Markov processes. For instance, in [11], one can find a case 
of the two-steps example (4 x 4). The second-order differential operator is also the infinitesimal 
operator of a bivariate Markov process. In this case, apart from the parameters a,(3 > —1, there 
are two new extra parameters, ki and /c2 with the condition that 0<A;i < k2 </3 + l. Therefore 
this example can be viewed as another variant of the Wright-Fisher model involving only mutation 
effects, but with two parameters to study, which decide how the process move through all the 
phases. 

Eigenfunctions of second-order differential operators with matrix- valued coefficients have been 
studied in the context of matrix- valued orthogonal polynomials, see [5l [TOl 111] . These examples 
are potential candidates of bivariate Markov processes. The difference is that these differential 
operators do not exactly fit into the framework of bivariate Markov processes. An inverse transfor- 



mation of the type (4.7) is going to be needed in order to get an infinitesimal operator of a bivariate 



Markov process, which in principle is not an easy task. This consideration, however, goes beyond 
the scope of this paper, and will be pursued elsewhere. 
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